---
title: "Table 4 results: Pancreatectomy"
params:
  procedure: pancreat
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE, cache = FALSE}
require("knitr")
## setting working directory for output
opts_knit$set(root.dir = "C:/Users/kcha0642/Documents/NSW-ON-NY/on-ny-nsw-kcha0642/on-ny-nsw")

opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, cache=TRUE, 
               comment=NA, verbose=TRUE, fig.width=6, fig.height=4, echo = FALSE)
```

```{r}
source("table4/table4 functions.R")
library(haven)
```


```{r}
# read cohort in
cohort <- read_sas(paste0(data_directory,"/",params$procedure,"cohort.sas7bdat"),
                   col_select = c(ppn,recnum,age,av_proc_year,SA2_2011_CODE,sex,
                                  los_eoc,index_date,death_date,
                                  mode_of_separation_recode,readmit_30,end_eoc,start_eoc,
                                  readmit_90,episode_start_date,episode_end_date,
                                  death_episode,death_episode_7days,discharge_home,
                                  discharge_longterm,discharge_other))
cohort_work <- cohort
```

```{r}
## table 1 check
print(sum(cohort$episode_start_date != cohort$start_eoc | cohort$episode_end_date != cohort$end_eoc))
print(sum(cohort$episode_start_date != cohort$start_eoc | cohort$episode_end_date != cohort$end_eoc)/nrow(cohort)*100)
```

```{r}
# set up data / transformations
# adjustments for ON mean values
cohort_work <- switch(params$procedure,
                 "pancreat" = adjust_covariates(cohort_work,mean.ontario.age = 64.81,
                                                mean.ontario.volume = 32.1),
                 "prostat" = adjust_covariates(cohort_work,mean.ontario.age = 62.76,
                                               mean.ontario.volume = 47.2),
                 "neph" = adjust_covariates(cohort_work,mean.ontario.age = 61.71,
                                            mean.ontario.volume = 39.4)
                 )
cohort_work <- sa2_income(cohort_work,sa2_cor)
```

```{r}
cohort_work <- binary_outcomes(cohort_work)
```

### Unadjusted outcomes
```{r}
los <- los_unadjusted(cohort_work)
print(los)
```

```{r}
binary_variables = c("death_episode","death_episode_7days","discharge_home","discharge_longterm",
                     "discharge_other","readmit_30","readmit_90")
n_binary = length(binary_variables)

bin_out = lapply(binary_variables,function(i){binary_unadjusted(cohort_work,i)})

for (i in 1:n_binary){
  print(binary_variables[i])
  print(bin_out[[i]])
}
```

### Adjusted outcomes
```{r}
cohort_work <- cohort_work %>%
  mutate(ppn = as.factor(ppn),
         sex = as.factor(sex)) %>%
  mutate_at(vars(binary_variables),as.numeric)
```


#### Model 1

##### length of stay
```{r}
M1_los <- model1_los(cohort_work,params$procedure)
print(summary(M1_los))
```

##### binary variables
```{r}
M1_bin <- lapply(binary_variables,function(i){model1_bin(cohort_work,i,params$procedure)})

for (i in 1:n_binary){
  print(binary_variables[i])
  print(M1_bin[[i]])
}
```

#### Model 2
##### length of stay
```{r}
M2_los <- model2_los(cohort_work,procedure=params$procedure)
print(summary(M2_los))
```

##### binary variables
```{r}
M2_bin <- lapply(binary_variables,function(i){model2_bin(cohort_work,i,params$procedure)})

for (i in 1:n_binary){
  print(binary_variables[i])
  print(M2_bin[[i]])
}
```

# Table outputs
```{r}
table4 <- data.frame(variable = c("LOS","LOS median","LOS prop","Died in hospital","Died within 7 days",
                                  "Home","Longterm","Other","Readmission 30","Readmission 90"),
                     unadjusted = NA,
                     model1 = NA,
                     model2 = NA)

table4$variable <- as.character(table4$variable)

mean_sd <- function(meanX,sdX){paste0(round(meanX,1)," (",round(sdX,1),")")}

## unadjusted values
table4 <- table4 %>%
  mutate(unadjusted = case_when(variable=="LOS" ~ mean_sd(los$los_mean,los$los_sd),
                                variable=="LOS median" ~ paste0(los$los_median," (",los$los_iqr,")"),
                                variable=="LOS prop" ~ as.character(los$los_100perc),
                                variable=="Died in hospital" ~ mean_sd(bin_out[[1]]$var.sum,
                                                                         bin_out[[1]]$var.prop*100),
                                variable=="Died within 7 days" ~ mean_sd(bin_out[[2]]$var.sum,
                                                                       bin_out[[2]]$var.prop*100),
                                variable=="Home" ~ mean_sd(bin_out[[3]]$var.sum,
                                                           bin_out[[3]]$var.prop*100),
                                variable=="Longterm" ~ mean_sd(bin_out[[4]]$var.sum,
                                                               bin_out[[4]]$var.prop*100),
                                variable == "Other" ~ mean_sd(bin_out[[5]]$var.sum,
                                                                          bin_out[[5]]$var.prop*100),
                                variable == "Readmission 30" ~ mean_sd(bin_out[[6]]$var.sum,
                                                                       bin_out[[6]]$var.prop*100),
                                variable == "Readmission 90" ~ mean_sd(bin_out[[7]]$var.sum,
                                                                       bin_out[[7]]$var.prop*100))  )

## adjusted model 1 values
adj_ci <- function(mod){
  m <- summary(mod)$coefficients["(Intercept)","Estimate"]
  se <- summary(mod)$coefficients["(Intercept)","Std.err"]
  return(paste0(round(m,1)," (",round(m-qnorm(0.975)*se,1),", ",round(m+qnorm(0.975)*se,1),")"))
}

adj_ci_b <- function(mod){
  if (is.character(mod)){
    val = ""
  } else{
    m <- summary(mod)$coefficients["(Intercept)","Estimate"]
    ci <- m + c(-1,1)*summary(mod)$coefficients["(Intercept)","Std.err"]*qnorm(0.975)
    m <- round(exp(m)/(1+exp(m))*100,1)
    ci <- round(exp(ci)/(1+exp(ci))*100,1)
    val <- paste0(m," (",ci[1],", ",ci[2],")")
  }
  return(val)
}

table4 <- table4 %>%
  mutate(model1 = case_when(variable=="LOS" ~ adj_ci(M1_los),
                                variable=="Died in hospital" ~ adj_ci_b(M1_bin[[1]]),
                                variable=="Died within 7 days" ~ adj_ci_b(M1_bin[[2]]),
                                variable=="Home" ~ adj_ci_b(M1_bin[[3]]),
                                variable=="Longterm" ~ adj_ci_b(M1_bin[[4]]),
                                variable == "Other" ~ adj_ci_b(M1_bin[[5]]),
                                variable == "Readmission 30" ~ adj_ci_b(M1_bin[[6]]),
                                variable == "Readmission 90" ~ adj_ci_b(M1_bin[[7]])))

## adjusted model 2 values
table4 <- table4 %>%
  mutate(model2 = case_when(variable=="LOS" ~ adj_ci(M2_los),
                                variable=="Died in hospital" ~ adj_ci_b(M2_bin[[1]]),
                                variable=="Died within 7 days" ~ adj_ci_b(M2_bin[[2]]),
                                variable=="Home" ~ adj_ci_b(M2_bin[[3]]),
                                variable=="Longterm"  ~ adj_ci_b(M2_bin[[4]]),
                                variable == "Other" ~ adj_ci_b(M2_bin[[5]]),
                                variable == "Readmission 30" ~ adj_ci_b(M2_bin[[6]]),
                                variable == "Readmission 90" ~ adj_ci_b(M2_bin[[7]])))

write_csv(table4,paste0("table4_",params$procedure,".csv"),na = "")


```

```{r}
## Get output for Vicki's analysis
output <- read_csv("Copy of SEs for Vicki.csv")
output %>%
  mutate(Estimate = NA,
         Std.err = NA) %>%
  mutate_at(vars(Estimate,Std.err),as.numeric) %>%
  filter(procedure == "PD") %>%
  select(-procedure) ->output

fill_table_m1 <- function(var_outcome,mod){
  output <- output %>%
    mutate(Estimate = ifelse(outcome == var_outcome & !grepl("income",model), tryCatch(summary(mod)$coefficients["(Intercept)","Estimate"],error = function(e){NA}),Estimate)) %>%
    mutate(Std.err = ifelse(outcome == var_outcome &  !grepl("income",model), tryCatch(summary(mod)$coefficients["(Intercept)","Std.err"],error = function(e){NA}),Std.err))
  return(output)
}

fill_table_m2 <- function(var_outcome,mod){
  output <- output %>%
    mutate(Estimate = ifelse(outcome == var_outcome & grepl("income",model), tryCatch(summary(mod)$coefficients["(Intercept)","Estimate"],error = function(e){NA}),Estimate)) %>%
    mutate(Std.err = ifelse(outcome == var_outcome & grepl("income",model), tryCatch(summary(mod)$coefficients["(Intercept)","Std.err"],error = function(e){NA}),Std.err))
  return(output)
}

outcomes = unique(output$outcome)
n_outcomes = length(outcomes)

for (i in 1:n_outcomes){
  if (i == 1){
    output <- fill_table_m1(var_outcome = "LOS",mod = M1_los)
    output <- fill_table_m2(var_outcome = "LOS",mod = M2_los)
  } else {
    output <- fill_table_m1(var_outcome = outcomes[i],mod = M1_bin[[i-1]])
    output <- fill_table_m2(var_outcome = outcomes[i],mod = M2_bin[[i-1]])
  }
}

output <- output %>%
  mutate(procedure = params$procedure)

output %>%
  write_csv("SEs for Vicki NSW.csv",append = TRUE,na = "")

```















